%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 8 Panel A: Leverage effect and the variance risk premium
% It also requires 1 function file: olshac.m
% Default code generates the case of correlation rho=0
% For the other correlation change rho to 0.415 in row 21
% Typically it takes less than 3 minutes to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;   % vary correlation to rho=0.415
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);

vec_corr=zeros(1,N_th);

for i=1:N_th
    theta=vec_theta(i);    % Vary theta
    % Implied Quantities after fixing theta
    M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
    M_c=theta/(alpha*(1+theta));
    M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
    M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
    M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
    M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
    M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
    M_corr_sv=M_Cov_con/sqrt(M_Var_con);
    vec_corr(i)=M_corr_sv;    %   Correlation
end
Correlation=round(vec_corr,2) % Report correlation (leverage effect)


disp('Press any key to continue to Predictive regressions ...');
pause;

% Simulations and regressions
h=3;       % predictive regression holding return horizon in months 1,3,12
N_sim=1000;                          % number of simulations: use 1000
S_end_year=1000;                     % sample size
int = 1/252;           %             % time interval for each year in a path
time_grid = 0:int:S_end_year;        % Time vector from 0 to end year
N_time=length(time_grid);
N_rets=N_time-1;                    % Number of returns in each path
H=1/12;     % Data frequency in the regression e.g. 1/12 for monthly
N_pp=H/int; % number of data points for each data point in regression: 21 for month


% omitting the first burn out period
keep_final_years=S_end_year; % keep all of it
Obs=keep_final_years*12;

% 1 path vectors
vec_D = zeros(1,N_time); % Allocate output vector
vec_V = zeros(1,N_time); % Allocate output vector
vec_X = zeros(1,N_time); % Allocate output vector

% All paths monthly matrices
N_month=(12*S_end_year)+1;
SM_V = zeros(1,N_month); % Allocate output vector
SM_X = zeros(1,N_month); % Allocate output vector
SM_D = zeros(1,N_month); % Allocate output vector
SM_S = zeros(1,N_month); % Allocate output vector
SM_VSR = zeros(1,N_month); % Allocate output vector
SM_VRP = zeros(1,N_month); % Allocate output vector
SM_RV  = zeros(1,N_month); % Allocate output vector

% Univariate regression of VRP
vec_beta_VRP        =zeros(N_sim,N_th);
vec_tstat_hac_VRP  =zeros(N_sim,N_th);
vec_R2_adj_VRP      =zeros(N_sim,N_th);

% Begin simulation and regression
rng('default'); % fix the random number generator  seed   for replication
vec_Z1=randn(N_sim,N_time);
vec_Z2=randn(N_sim,N_time);
tic
for s = 1:N_sim
    for th=1:N_th
        theta=vec_theta(th);
        % Stock price quantities
        M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
        M_c=theta/(alpha*(1+theta));
        M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
        M_a=a;
        M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
        M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
        M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
        M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
        
        
        % VSR and VRP quantities
        M_zeta_v=zeta*M_Var_con;
        M_kappa_v=kappa+(sigma*M_Cov_con);
        M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
        v0=M_zeta_v/kappa; % LR stock variance level
        M_RN_LRM_vt=(M_zeta_v+(r*M_eta_v))/(M_kappa_v+(0.5*M_eta_v));
        
        % Initalization of each path
        D0=100;
        V0=zeta/kappa; % long run mean of  fundamental  variance
        X0=(r+(0.5*M_Var_con*V0))/(1+theta); %  inital  X is equal to its long run mean
        
        vec_D(1)=D0;
        vec_V(1)=V0;
        vec_X(1)=X0;
        
        for j = 1:N_rets
            Z1=vec_Z1(s,j);
            Z2=vec_Z2(s,j);
            vec_V(j+1) = max(0, vec_V(j))+((zeta-(kappa*max(0, vec_V(j))))*int)...
                +(sigma*rho*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*int)...
                +(alpha*M_sigma_1*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1)...
                +(alpha*M_sigma_2*sqrt(max(0, vec_V(j)))*sqrt(int)*Z2);
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*int)+(vec_D(j)*sqrt(max(0, vec_V(j)))*sqrt(int)*Z1);
            
        end
        
        % Get end of month values for each process from daily values
        SM_V= vec_V(1:N_pp:end);
        SM_X= vec_X(1:N_pp:end);
        SM_D= vec_D(1:N_pp:end);
        SM_S= SM_D.*exp(M_a+(M_b*SM_V)+(M_c*SM_X)); %stock price
        
        % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################      instantenous  RV       ################
        %  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        RV_A= M_zeta_v;
        RV_B=1-kappa;
        
        % +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        %  #######################   instantenous     VSR      ################
        %  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        VSR_A=M_zeta_v;
        VSR_B=1-M_kappa_v;
        VSR_C=M_eta_v;

        % Simulated VSR VRP and RV
        SM_RV=(RV_A+(RV_B*M_Var_con.*SM_V))*(10000/12); %normalized to get monthly values in percentage squared
        SM_VSR=(VSR_A+(VSR_B*M_Var_con.*SM_V)+(VSR_C.*SM_X))*(10000/12);  %normalized to get monthly values in percentage squared
        SM_VRP=SM_RV-SM_VSR;        
        
        % Keep only the selected final years and omit the first years as burnout if needed
        SM_RV=SM_RV(end-Obs:end);
        SM_VRP=SM_VRP(end-Obs:end);
        SM_S=SM_S(end-Obs:end);
        SM_D=SM_D(end-Obs:end);
        
        % Get h period ahead log returns
        SM_Ret_h=log(SM_S(1+h:end))-log(SM_S(1:end-h)); % log ret
        %SM_Ret_h=(SM_S(1+h:end)./SM_S(1:end-h))-1;      % returns
        SM_Ret_h_ann=(12/h)*SM_Ret_h*100;  %annualized in percentage terms
        
        
        % Univariate Regression of h period returns on VRP
        % ===============================================
        
        reg_Y=SM_Ret_h_ann';            %   annualized   returns in %
        LT=length(reg_Y);
        reg_X=[ones(LT,1) SM_VRP(1:end-h)'];
        lag=round(0.75*(LT^(1/3)));    %  recommended hac lag choice
        
        [beta, R2, R2adj, X2_NW, X2_HH, X2_R, std_NW, std_HH, std_R, t_NW, t_HH, t_R]=olshac(reg_Y,reg_X,lag,lag);
        
        vec_beta_VRP(s,th)=beta(2);
        vec_tstat_hac_VRP(s,th)=t_NW(2);
        vec_R2_adj_VRP(s,th) = R2adj;       
        
    end   % end theta
    
end    % end sim
clc;
toc

%  REPORT   UNIVARIATE REGRESSION on VRP QUANTITIES
tstat_hac_VRP=round(mean(vec_tstat_hac_VRP),2)
R2_VRP=round(mean(vec_R2_adj_VRP)*100,2)


clear alpha D0 EE0 EE1 EE2 EEv j jj kappa kappa1 kappa2 m M_a M_b M_c;
clear M_corr_sv M_Cov_con M_DELTA M_eta_v M_kappa_v M_Ob_LRM_vt M_r M_RN_LRM_vt;
clear M_sigma_1 M_sigma_2 M_sigma_v  M_Var_con M_zeta_v mu n rho sigma;
clear theta tt1 v V0 Z1 Z2 zeta;
clear vec_Z1 vec_Z2;

%% THE END










